/* ssyevx.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
 */

#include "Lapack.h"
#include <cmath>

using namespace std;

extern "C" {
	int lsame_(char*, char*);
	int xerbla_(char *, int *);
}

extern double slamch_(char *);
extern int ilaenv_(int *, const char *, const char *, int *, int *, int *, int *);

static int c__1 = 1;
static int c_n1 = -1;

int Lapack::ssyevx(char *jobz, char *range, char *uplo, int *n,
		float *a, int *lda, float *vl, float *vu, int *il, int *iu,
		float *abstol, int *m, float *w, float *z__, int *ldz, float *
		work, int *lwork, int *iwork, int *ifail, int *info) 
{
	/* System generated locals */
	int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
	float r__1, r__2;

	/* Local variables */
	int i__, j, nb, jj;
	float eps, vll, vuu, tmp1;
	int indd, inde;
	float anrm;
	int imax;
	float rmin, rmax;
	int test;
	int itmp1, indee;
	float sigma;
	int iinfo;
	char order[1];
	int lower;
	int wantz, allEigen, indexEigen;
	int iscale, indibl;
	int valeig;

	float safmin;
	float abstll, bignum;
	int indtau, indisp, indiwo, indwkn;
	int indwrk, lwkmin;
	int llwrkn, llwork, nsplit;
	float smlnum;
	int lwkopt;
	int lquery;

	/*  -- LAPACK driver routine (version 3.2) -- */
	/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
	/*     November 2006 */

	/*     .. Scalar Arguments .. */
	/*     .. */
	/*     .. Array Arguments .. */
	/*     .. */

	/*  Purpose */
	/*  ======= */

	/*  SSYEVX computes selected eigenvalues and, optionally, eigenvectors */
	/*  of a float symmetric matrix A.  Eigenvalues and eigenvectors can be */
	/*  selected by specifying either a range of values or a range of indices */
	/*  for the desired eigenvalues. */

	/*  Arguments */
	/*  ========= */

	/*  JOBZ    (input) CHARACTER*1 */
	/*          = 'N':  Compute eigenvalues only; */
	/*          = 'V':  Compute eigenvalues and eigenvectors. */

	/*  RANGE   (input) CHARACTER*1 */
	/*          = 'A': all eigenvalues will be found. */
	/*          = 'V': all eigenvalues in the half-open interval (VL,VU] */
	/*                 will be found. */
	/*          = 'I': the IL-th through IU-th eigenvalues will be found. */

	/*  UPLO    (input) CHARACTER*1 */
	/*          = 'U':  Upper triangle of A is stored; */
	/*          = 'L':  Lower triangle of A is stored. */

	/*  N       (input) INTEGER */
	/*          The order of the matrix A.  N >= 0. */

	/*  A       (input/output) REAL array, dimension (LDA, N) */
	/*          On entry, the symmetric matrix A.  If UPLO = 'U', the */
	/*          leading N-by-N upper triangular part of A contains the */
	/*          upper triangular part of the matrix A.  If UPLO = 'L', */
	/*          the leading N-by-N lower triangular part of A contains */
	/*          the lower triangular part of the matrix A. */
	/*          On exit, the lower triangle (if UPLO='L') or the upper */
	/*          triangle (if UPLO='U') of A, including the diagonal, is */
	/*          destroyed. */

	/*  LDA     (input) INTEGER */
	/*          The leading dimension of the array A.  LDA >= max(1,N). */

	/*  VL      (input) REAL */
	/*  VU      (input) REAL */
	/*          If RANGE='V', the lower and upper bounds of the interval to */
	/*          be searched for eigenvalues. VL < VU. */
	/*          Not referenced if RANGE = 'A' or 'I'. */

	/*  IL      (input) INTEGER */
	/*  IU      (input) INTEGER */
	/*          If RANGE='I', the indices (in ascending order) of the */
	/*          smallest and largest eigenvalues to be returned. */
	/*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
	/*          Not referenced if RANGE = 'A' or 'V'. */

	/*  ABSTOL  (input) REAL */
	/*          The absolute error tolerance for the eigenvalues. */
	/*          An approximate eigenvalue is accepted as converged */
	/*          when it is determined to lie in an interval [a,b] */
	/*          of width less than or equal to */

	/*                  ABSTOL + EPS *   max( |a|,|b| ) , */

	/*          where EPS is the machine precision.  If ABSTOL is less than */
	/*          or equal to zero, then  EPS*|T|  will be used in its place, */
	/*          where |T| is the 1-norm of the tridiagonal matrix obtained */
	/*          by reducing A to tridiagonal form. */

	/*          Eigenvalues will be computed most accurately when ABSTOL is */
	/*          set to twice the underflow threshold 2*SLAMCH('S'), not zero. */
	/*          If this routine returns with INFO>0, indicating that some */
	/*          eigenvectors did not converge, try setting ABSTOL to */
	/*          2*SLAMCH('S'). */

	/*          See "Computing Small Singular Values of Bidiagonal Matrices */
	/*          with Guaranteed High Relative Accuracy," by Demmel and */
	/*          Kahan, LAPACK Working Note #3. */

	/*  M       (output) INTEGER */
	/*          The total number of eigenvalues found.  0 <= M <= N. */
	/*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */

	/*  W       (output) REAL array, dimension (N) */
	/*          On normal exit, the first M elements contain the selected */
	/*          eigenvalues in ascending order. */

	/*  Z       (output) REAL array, dimension (LDZ, max(1,M)) */
	/*          If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
	/*          contain the orthonormal eigenvectors of the matrix A */
	/*          corresponding to the selected eigenvalues, with the i-th */
	/*          column of Z holding the eigenvector associated with W(i). */
	/*          If an eigenvector fails to converge, then that column of Z */
	/*          contains the latest approximation to the eigenvector, and the */
	/*          index of the eigenvector is returned in IFAIL. */
	/*          If JOBZ = 'N', then Z is not referenced. */
	/*          Note: the user must ensure that at least max(1,M) columns are */
	/*          supplied in the array Z; if RANGE = 'V', the exact value of M */
	/*          is not known in advance and an upper bound must be used. */

	/*  LDZ     (input) INTEGER */
	/*          The leading dimension of the array Z.  LDZ >= 1, and if */
	/*          JOBZ = 'V', LDZ >= max(1,N). */

	/*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
	/*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */

	/*  LWORK   (input) INTEGER */
	/*          The length of the array WORK.  LWORK >= 1, when N <= 1; */
	/*          otherwise 8*N. */
	/*          For optimal efficiency, LWORK >= (NB+3)*N, */
	/*          where NB is the max of the blocksize for SSYTRD and SORMTR */
	/*          returned by ILAENV. */

	/*          If LWORK = -1, then a workspace query is assumed; the routine */
	/*          only calculates the optimal size of the WORK array, returns */
	/*          this value as the first entry of the WORK array, and no error */
	/*          message related to LWORK is issued by XERBLA. */

	/*  IWORK   (workspace) INTEGER array, dimension (5*N) */

	/*  IFAIL   (output) INTEGER array, dimension (N) */
	/*          If JOBZ = 'V', then if INFO = 0, the first M elements of */
	/*          IFAIL are zero.  If INFO > 0, then IFAIL contains the */
	/*          indices of the eigenvectors that failed to converge. */
	/*          If JOBZ = 'N', then IFAIL is not referenced. */

	/*  INFO    (output) INTEGER */
	/*          = 0:  successful exit */
	/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
	/*          > 0:  if INFO = i, then i eigenvectors failed to converge. */
	/*                Their indices are stored in array IFAIL. */

	/* ===================================================================== */

	/*     .. Parameters .. */
	/*     .. */
	/*     .. Local Scalars .. */
	/*     .. */
	/*     .. External Functions .. */
	/*     .. */
	/*     .. External Subroutines .. */
	/*     .. */
	/*     .. Intrinsic Functions .. */
	/*     .. */
	/*     .. Executable Statements .. */

	/*     Test the input parameters. */

	/* Parameter adjustments */
	a_dim1 = *lda;
	a_offset = 1 + a_dim1;
	a -= a_offset;
	--w;
	z_dim1 = *ldz;
	z_offset = 1 + z_dim1;
	z__ -= z_offset;
	--work;
	--iwork;
	--ifail;

	/* Function Body */
	lower = lsame_(uplo, "L");
	wantz = lsame_(jobz, "V");
	allEigen = lsame_(range, "A");
	valeig = lsame_(range, "V");
	indexEigen = lsame_(range, "I");
	lquery = *lwork == -1;

	*info = 0;
	if (!(wantz || lsame_(jobz, "N"))) {
		*info = -1;
	} else if (!(allEigen || valeig || indexEigen)) {
		*info = -2;
	} else if (!(lower || lsame_(uplo, "U"))) {
		*info = -3;
	} else if (*n < 0) {
		*info = -4;
	} else if (*lda < max(1, *n)) {
		*info = -6;
	} else {
		if (valeig) {
			if (*n > 0 && *vu <= *vl) {
				*info = -8;
			}
		} else if (indexEigen) {
			if (*il < 1 || *il > max(1, *n)) {
				*info = -9;
			} else if (*iu < min(*n, *il) || *iu > *n) {
				*info = -10;
			}
		}
	}
	if (*info == 0) {
		if (*ldz < 1 || wantz && *ldz < *n) {
			*info = -15;
		}
	}

	if (*info == 0) {
		if (*n <= 1) {
			lwkmin = 1;
			work[1] = (float) lwkmin;
		} else {
			lwkmin = *n << 3;
			nb = ilaenv_(&c__1, "SSYTRD", uplo, n, &c_n1, &c_n1, &c_n1);
			/* Computing MAX */
			i__1 = nb, i__2 = ilaenv_(&c__1, "SORMTR", uplo, n, &c_n1, &c_n1,
					&c_n1);
			nb = max(i__1, i__2);
			/* Computing MAX */
			i__1 = lwkmin, i__2 = (nb + 3) * *n;
			lwkopt = max(i__1, i__2);
			work[1] = (float) lwkopt;
		}

		if (*lwork < lwkmin && !lquery) {
			*info = -17;
		}
	}

	if (*info != 0) {
		i__1 = -(*info);
		xerbla_("SSYEVX", &i__1);
		return 0;
	} else if (lquery) {
		return 0;
	}

	/*     Quick return if possible */

	*m = 0;
	if (*n == 0) {
		return 0;
	}

	if (*n == 1) {
		if (allEigen || indexEigen) {
			*m = 1;
			w[1] = a[a_dim1 + 1];
		} else {
			if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
				*m = 1;
				w[1] = a[a_dim1 + 1];
			}
		}
		if (wantz) {
			z__[z_dim1 + 1] = 1.f;
		}
		return 0;
	}

	/*     Get machine constants. */

	safmin = slamch_("Safe minimum");
	eps = slamch_("Precision");
	smlnum = safmin / eps;
	bignum = 1.f / smlnum;
	rmin = sqrt(smlnum);
	/* Computing MIN */
	r__1 = sqrt(bignum), r__2 = 1.f / sqrt(sqrt(safmin));
	rmax = min(r__1, r__2);

	/*     Scale matrix to allowable range, if necessary. */

	iscale = 0;
	abstll = *abstol;
	if (valeig) {
		vll = *vl;
		vuu = *vu;
	}
	anrm = slansy("M", uplo, n, &a[a_offset], lda, &work[1]);
	if (anrm > 0.f && anrm < rmin) {
		iscale = 1;
		sigma = rmin / anrm;
	} else if (anrm > rmax) {
		iscale = 1;
		sigma = rmax / anrm;
	}
	if (iscale == 1) {
		if (lower) {
			i__1 = *n;
			for (j = 1; j <= i__1; ++j) {
				i__2 = *n - j + 1;
				sscal(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
			}
		} else {
			i__1 = *n;
			for (j = 1; j <= i__1; ++j) {
				sscal(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
			}
		}
		if (*abstol > 0.f) {
			abstll = *abstol * sigma;
		}
		if (valeig) {
			vll = *vl * sigma;
			vuu = *vu * sigma;
		}
	}

	/*     Call SSYTRD to reduce symmetric matrix to tridiagonal form. */

	indtau = 1;
	inde = indtau + *n;
	indd = inde + *n;
	indwrk = indd + *n;
	llwork = *lwork - indwrk + 1;
	ssytrd(uplo, n, &a[a_offset], lda, &work[indd], &work[inde], &work[
			indtau], &work[indwrk], &llwork, &iinfo);

	/*     If all eigenvalues are desired and ABSTOL is less than or equal to */
	/*     zero, then call SSTERF or SORGTR and SSTEQR.  If this fails for */
	/*     some eigenvalue, then try SSTEBZ. */

	test = false;
	if (indexEigen) {
		if (*il == 1 && *iu == *n) {
			test = true;
		}
	}
	
	if ((allEigen || test) && *abstol <= 0.f) {
		scopy(n, &work[indd], &c__1, &w[1], &c__1);
		indee = indwrk + (*n << 1);
		if (!wantz) {
			i__1 = *n - 1;
			scopy(&i__1, &work[inde], &c__1, &work[indee], &c__1);
			ssterf(n, &w[1], &work[indee], info);
		} else {
			slacpy("A", n, n, &a[a_offset], lda, &z__[z_offset], ldz);
			sorgtr(uplo, n, &z__[z_offset], ldz, &work[indtau], &work[indwrk]
					, &llwork, &iinfo);
			i__1 = *n - 1;
			scopy(&i__1, &work[inde], &c__1, &work[indee], &c__1);
			ssteqr(jobz, n, &w[1], &work[indee], &z__[z_offset], ldz, &work[indwrk], info);
			if (*info == 0) {
				for (i__ = 1; i__ <= *n; ++i__) {
					ifail[i__] = 0;
				}
			}
		}
		if (*info == 0) {
			*m = *n;
			goto L40;
		}
		*info = 0;
	}	

	/*     Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN. */

	if (wantz) {
		*(unsigned char *) order = 'B';
	} else {
		*(unsigned char *) order = 'E';
	}
	indibl = 1;
	indisp = indibl + *n;
	indiwo = indisp + *n;
	*m = *n;
	sstebz(range, order, n, &vll, &vuu, il, iu, &abstll, &work[indd], &work[
			inde], m, &nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[
			indwrk], &iwork[indiwo], info);

	if (wantz) {
		sstein(n, &work[indd], &work[inde], m, &w[1], &iwork[indibl], &iwork[
				indisp], &z__[z_offset], ldz, &work[indwrk], &iwork[indiwo], &
				ifail[1], info);

		/*        Apply orthogonal matrix used in reduction to tridiagonal */
		/*        form to eigenvectors returned by SSTEIN. */

		indwkn = inde;
		llwrkn = *lwork - indwkn + 1;
		sormtr("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau], &z__[
				z_offset], ldz, &work[indwkn], &llwrkn, &iinfo);
	}

	/*     If matrix was scaled, then rescale eigenvalues appropriately. */

L40:
	if (iscale == 1) {
		if (*info == 0) {
			imax = *m;
		} else {
			imax = *info - 1;
		}
		r__1 = 1.f / sigma;
		sscal(&imax, &r__1, &w[1], &c__1);
	}

	/*     If eigenvalues are not in order, then sort them, along with */
	/*     eigenvectors. */

	if (wantz) {
		i__1 = *m - 1;
		for (j = 1; j <= i__1; ++j) {
			i__ = 0;
			tmp1 = w[j];
			i__2 = *m;
			for (jj = j + 1; jj <= i__2; ++jj) {
				if (w[jj] < tmp1) {
					i__ = jj;
					tmp1 = w[jj];
				}
			}

			if (i__ != 0) {
				itmp1 = iwork[indibl + i__ - 1];
				w[i__] = w[j];
				iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
				w[j] = tmp1;
				iwork[indibl + j - 1] = itmp1;
				sswap(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1], &c__1);
				if (*info != 0) {
					itmp1 = ifail[i__];
					ifail[i__] = ifail[j];
					ifail[j] = itmp1;
				}
			}
		}
	}

	/*     Set WORK(1) to optimal workspace size. */

	work[1] = (float) lwkopt;

	return 0;

	/*     End of SSYEVX */

} /* ssyevx_ */
